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1 . Introduction 

Linear matrix equations play n very important role in system theory. In 

this paper we undertake the study of linear matrix equations which take the form 
s t 

Q (1.1) 

i-0 j»0 


I X “i/ p aJ ' 


where B(mxm), A(nxn) and Q(mxn) are given matrices over some field F and 
gjj are elements of F, using methods of modern algebra. The emphasis is on 
the solution of such equations using finite algebraic procedures which are 
easily implemented on a digital computer. 

Particular attention is given to equations PA + BP ■ Q and P - BPA ■ Q, 
and their special subcases PA + A'P ■» Q (the Lyapunov equation) and P - A' PA - Q 
the (discrete Lyapunov equation). The Lyapunov equation appears in several areas 
of control theory such ns stability theory, optimal control (evaluation of quad- 
ratic integrals), stochastic control (evaluation of covariance matrices) and in 
the solution of the Algebraic Riccati Equation using Newton's Method. 

This paper has been inspired by an important paper by Kalman [2]. Kal- 
man's concern was the characterization of polynomials whose zeros lie in 
certain algebraic domains (and the unification of the ideas of Hermite and 
Lyapunov). In this paper we show that the same ideas lead to finite algorithms 
for the solution of linear matrix equations of the form given above. The Analysis 
in terms of a module structure on matrices presented here is believed to be new. 

In a subsequent paper we shall investigate the implications of these ideas on 
stability theory. 

The paper is divided into five sections. In section 2 we define the action 


fjj^ over an arbitrary commutative ring with identify and prove a Basic 
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Lemma. In section 3 we consider equation (1.1) over a field F in great 
generality and prove the Main Theorem. In section 4 wo deal with the 
equation PA + BP ■ Q and the Lyapunov Equation PA + A*P « Q for which we 
give algorithms for obtaining its solution and comment on the arithmetic 
complexity. We also provide numerical examples and prove a stability 
theorem. In section 5 we deal with the equation P - BPA ■ Q as well as 
with the Discrete Lyapunov Equation P - A' PA « y. In section 6 we look 
at equation (1.1) over an integral domain. 

2. The Action f . 

BA 

Let A be an nxn matrix and B an mxm matrix both over E, a commutative 
ring with identity. Let E[x,y] be the ring of polynomials in two indetcr- 
minatos x and y over E. Let 4* = (<J> 2 (x), ^(y)) be t.he ideal in E[x,y] 
generated by $ 2 (x) the characteristic polynomial of A, and ^(y) the 
characteristic polynomial of B. Elements of the quotient ring E[x,y]/V 
are cosets (equivalence classes) denoted by Y + a(x,y). The Cayley 
Hamilton Theorem holds (41 therefore <f> 2 (A) = 0,^(0) = 0. 

Since <> 2 (x) and <^ 2 (y) ar *' monic polynomials division is possible 
and as a consequence we can state ; 

L emma 1 ; Let g(x,y)e E(x,yl. Then g(x,y) can be written uniquely ass 
g(x,y) = t (x,y) <J> 2 (x) ^(y) + p(x,y)<f> 2 (x) + q(x,y)ij> 2 (y) + r(x f y) 
where : 

the degree of p(x,y) in y is less than m (it may be a polynomial in x) or 


p(x,y) is zero, 
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tho degree of q(x,y) in x ie less than n (it may bo a (2.1) 

polynomial in y) or q(x,y) is zero, 

the degree of r(x,y) in y is loss than m, in x less than 
n or r(x,y) is zero- 

Proof t 

Division in x by ^ (x) is possible therefore 
g(x,y) - a(x,y)4> 2 (x) + b(x,y) 

where 

degree of b(x,y) in x is less than n (b(x,y) may be a polynomial in y) 
or b(x,y) is zero. 

Division in y by ^ 2 (y) is possible therefore 
a (x,y) = t(x,y)^ 2 (y) + p(x,y) 

where degree of p(x,y) in y is less than m (p(x,y) may be a polynomial 
in x) or p(x,y) is zero. 

Also 

b(x,y) = q(x,y)i|^ 2 (y) + r(x,y) 

where degree of r(x,y) in y is less than m and degree of r(x,y) in x is 
less than n or r(x,y) is zero. 

Now then 

g(x,y) = t (x,y)c}> 2 (x) i^ 2 (y) + p(x,y)<}> 2 (x) + q(x,y)v^ 2 (y) + r(x,y) 


This representation is unique since suppose 


g(x,y) 


- t^ (x,y) <f> 2 (x) <|< 2 (y) + P x (x,y) <> 2 vx) + q x (x,y) ^ (y) +r 1 <x,y) 

■ t 2 (x, y) ^ (x) <^ 2 (y) + P 2 (x,y) $ 2 (x) + q 2 (x,y) ^ (y) + r 2 (x,y) 

with p, , p , q , q , r , r satisfying requirements (2.1). 

r x (x,y) - r 2 (x,y) - [t^tj] ♦j(x) ^(y) + IPj-PjHjtx) + tq i -q 2 J ^ 2 (y) 

a BY 

Suppose that 01 ■/ 0 . Then there exists a term on the r.h.s. say a x* 

i > n, j >n. This term cannot be cancelled by either P or Y. Therefore, 
a ■ 0. Suppose that B ^ 0. Then there exists a term on the r.h.s. say 

b x* y^ i > n. This term cannot be cancelled by any term from Y. There- 

fore B ■ 0. But then 6 = 0 as well and r^(x,y) - r 2 (x,y). g 


Corollary 1 : Let g x - tj 4> 2 (x) ^ (y) + p^U) + q^ty) + and g 2 - 

t 2 <> 2 (x) ^(y) + P 2 4> 2 (x) + q 2 ^ 2 (y) + r 2 be in the same coset V + a(x,y). 
Then r, ■ r_. (If g - t<f> + p<P. + q4»_ + * denote r by g(x,y)mod f.) 

12 Z 2 2 2 

Let MW be the set of mxn matrices over E. Define the action 

f : E(x,y]x MW -*• MW in the following manner: 

BA 


f BA (h(X,y) ' M) 


where h(x,y) ■ £ 

jk 

MW. 




h B^ M A k 
3k 


jk 

is an element in Elx,yJ 


and M an element in 


It can be shown that f^ has the following properties. 

i) f (u,M) = uM where u e E. 

BA 

ii) f BA ^(x,y)+h(x,y) , M) = f Bft (g(x,y) ,M) + f^thtx.y) ,M) 
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iii) f BA (g(x,y)h(x,y) , M) - f BA (g(*»y)* f BA * h * x » y *' M) > 

- f BA ( h ( x »y) » f B A tg(x ' y) ' M)) 

iv) f BA (g( x *y)» M > " f BA (g(x,y)mod % M) 

v ) t BA (g(x ' y) ' M+N) “ f BA (g(X ' y) * M) + f BA (g(X ' y) ' N) 

Properties i) , ii), iii) and v) follow directly from the definition 

of f [1] . Property v) is arrived at by using Lemma 1 and the Cayloy- 
BA 

Hamilton theorem. 

..yt 

The definition of f _ allows for the interpretation of MN as an 

BA 

E{x,y]/V-module. 

Basic Lemma : The set MM of mxn matrices with elements in E is a module 

over the quotient ring Elx,y ]/¥• 

Proof : The set of mxn matrices under addition is an abelian group. Define 

multiplication (*) of cosets ¥ + h(x,y) and mxn matrices M in the follow- 
ing manners 

(’i'fh (x,y» * M « f BA (h(x,y)mod VjM) 

The multiplication is well defined and satisfies the properties: 

1) (4'+h(x,y)) *(M+N) = (y+h (x,y)) *M + (Y+h (x,y) ) *N 

2) (Y+h (x,y)) *[ (4'+g (x,y) *M] » [( , J'+h(x,y) ) • (4 +g (x,y) ) ] *M 

3) ( (V+h (x,y) ) + (4 / +g (x,y) ) ] *M = W+Mx.y) )*M + (f+g(x,y))*M 

4) (T+D^M - M 

for all M,N in MM and all T+h(x,y), T+gfXjy) in ETx^lA' with Y+l being 
the multiplicative identity in Etx.yl/i'. 
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3. The Ccneral Equation 

Suppose that we restrict E to bo some fiold F and lot K bo an alge- 
braically closed extension of F. If f(x,y) is an element of F[x,y) we 

K (1) 

denote by V f the variety of f(x,y) in A* 17). Let X^^ X^,..., n be the 

eigenvalues of A and p, , the eigenvalues of D. Suppose that 

12 m 

* 1 k 

g(x,y) is a polynomial in Flx,y). If g(x,y)“ \ g^ k y x then we define 


G the mnxmn matrix. 
9 


% '1 g jk b3 © a ' 


(3.1) 


where (x} denotes tensor product, (A(^B - (n„ B) ) and A' denotes 
transpose. The significance of the matrix comes from the following. 

Let £ be the mnxl column vector made up of the entries of matrix 
P - (p„) written as 

E ■ 'Ell P 12 P 13 — P ln P 21 P 22— P 2n p n l P ™2 ' ' '^r, 1 ' ' 

Let £ be the mnxl column vector made up of the entries of Q. Then 
equation (1.1) can simply be written as 


G g £ - a 


(3.2) 


Wo now state the 


Main Theorem : The following statements are equivalent. 

1) Equation (1.1) has a unique solution for all Q. 

2) G is invertible. 

9 

3) g(X i ,Pj) / 0 VnA l<i<n l<j<m. 


- { (t x » t 2 ) | t lt t 2 e K). 


(1) A* 


4 ) 


v g(x,y) O v $ 2 (x) nv ij/ 2 (y) " * 

5) The cooet 'f' , +g(x,y) is a unit in F[x,yJA» 

Proof : 

Wo will show the equivalences in the order 1) 2) 3) ^ 4)-=^ 5)^ 1) 

1)^ 2 

Suppose then that equation (1.1) does have a unique solution for all Q. 

Well since equation (1.1) can equivalently be written as G p - q then G 

9 9 

is invertible. 

2) -^3) 

From (5] theorem 43.8 we have that the A characteristic values of 

n 

G^ are gJA^P.j). Since det |"g(A i# Vi.. ) and since dot G^ / 0 we have 

that g(V ,p_.) 1 0 for all A ^ , v» ^ 1 < i < n 1 < j < m. 

3 ) ^ 4 ) 

If we look at what 3) says it is the following; that the polynomials 
g (x,y) , <j> 2 (x) and ^ (/) have no comnan roots in A^. But this is statement 

4) . 

4)=^75) 

Now Y +g(x,y) is a unit iff there exists a 4'+f(x,y) such that 

(*tf(x,y)M'F+g(x,y)) - ^ + 1 . 

Now 

4 t tg(x,y) is a unit <-^> 3 f(x,y),3 V + f (x, y)g (c f y) “ Y + 1 

< > 3 f(x,y), a 1 (x,y), a 2 (x,y) c F[x,y] such that 

f(x,y)g(x,y) + a^ (x»y)<|> 2 (x) + a 2 (x,y)ij> 2 (y) - 1 
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Assume now that 4) holds (i.o., the polynomial h»l vanishes at every com- 
mon zero of g (x,y) , <> 2 (x) , >J/ 2 (y) .) By tho Hilbcrt-Nullstcllonsatz (71 there 
exist polynomials f(x,y), a^x^) a 2 (x,y) such that 

f(x,y) g(x,y) + a A (x,y)4* 2 (x) + a 2 (x,y) ^ (y) - 1 


this means that Y+g(x,y) is a unit in F[x,y]/Y. 

5) ^1) 

Suppose that Y+g(x,y) is a unit in F(x,y)/Y i.c.J Y+f (x,y) such that 

(V*f (x.y)MV+g(x,y))- Y + 1. Lot P - f BA < f <x,y)mod Y, Q) - t^itlx.y) . Q) . 

Show that this is a solution to (1.1). 
s t 

X X 9 ij fil P A> ■ f BA <9(X ' y) ' P) 

i“0 j-0 


- f BA (g(x,y)f (x,y), g) 

■ W l - S) - » 


The P so defined is the unique solution to (1.1). 

Let P^ f P 2 be two distinct solutions to (1.1) 

^ f BA (9(X ' y) ' P l ) “ f BA (9(X ' y) ' P 2 > “ C 

f BA (f(X ' y>r f BA <g(X ' y) ' P l )) = f BA <f(X ' y) ' t BA (g(X ' y) ' P 2 )J 
P 1 ’ P 2 

which is a contradiction. 


Therefore equation (1.1) has a unique solution for all Q. This completes 
the proof of the Main Theorem, era 
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Romark 1 In the above proof wo havo an explicit expression for tho 
solution of equation (1.1). A general method for constructing such an 
f(x,y) is through a constructivo proof of tho Ililbcrt-Nullstcllcnsatz or 
using Resultant Theory (C) . As will bo soen in lator pages of this paper 
for several important equations this generality is unnecessary and easier 
methods exist. 

R emark 2 In our entiro construction we have been using the ideal 

V - (4> 2 (x) , (y) ) • Other ideals can be used. As an example tho ideal 

($ 2 (x) , ^ (y) ) where 4 2 (x) an< ^ ^(y) aro the minimal polynomicals of A 

and B respectively. Since ^(x) ■ k(x)<^(x) and ^(y) » My)<P 2 (y) we 

will be dealing with polynomials of smaller degree. This may have as an 

effect the reduction in the number of computations performed. 

Remark 3 In the special case in which A is missing from equation (1.1) 

« i 

(i.e., suppose it is of the form g. B P ■ Q) then it would seem 

i=0 

that the analysis can take place in some quotient ring F[y]/4>. This 

actually is the case. Lot 4 be the ideal in F[yJ generated by ^(y). Then 

¥+g(y) is a unit in F(x,yl/f if and only if 'Hg(y) is a unit in F [y] /<!>. 

This follows from the fact that if 3 fly), a 2 (y) elements of F(y) such 

that f(y)g(y) + a 2 (y)i^ 2 (y) *» 1 then clearly there exist elements 

f (x,y) (=f (y) ) , a (x,y) (-0) and a 2 (x,y) (=a 2 (y) in Flx,y] such that 

f(x,y)g(y) + (x,y) <f> 2 (x) + a 2 (x,y) i|> 2 (y) ■ 1. On the other hand if there 

exist f(x,y), a^fx.y), a 2 (x,y) elements of F[x,y] such that f(x,y)g(y) + 

a^(x,y) x n + a 2 (x,y) \p 2 (y) ■ 1 then evaluating at x=0 we get f(o,y)g(y) + 

a 2 (0,y)4> 2 (y) ■ 1 which means that 4+g(y) is a unit in Fly] /4. (4> 2 (x) ■ 

x n since A = 0 ). The action f x F[y]x MM -*• MM can similarly be 
n B 
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defined 03 f (h (y) , M) h B^P and MW becomes ar. HyJ/ 4. -module. The 
B *"1 J 

® j J 

solution to'' g. B P * Q is then given by P ■ f n (f(x)mod Pi . We act 

I H 

i -0 

similarly if B is missing. 

Rema rk 4 Let U3 look at the vory special case wh_i» we are dealing with 

equation B£ ■ q whore £ and g[ are mxl vectors. In this case g(x,y) « y. 

What we want to do is find f (y) such that f(y)y + a(y)^ 2 (y) - 1. If 

♦_(y) - y m + k .y m_1 +...+k obvious choices are f(y)““-r^- y m 1 - 
2 m-1 o Xq 

Up y m ~ 2 - ••• - zr 1 a(y) “ k" 8ince 

K o 0 K o 

K )c v if 

• 1 m-1 m-1 m- 2 1 1 m m— 1 ro— 1 1 

<- r- y - r~ y - ••• - z^y + y + v y ♦ ... + ^ y + 1 - 1. 

0 K o K o k o k o k 0 

Now k / 0 since for a solution to exist detB - k ft / 0. The 


solution £ is given by: 


£ - f B (f(y), a) 


Analyzing this further wo get 

m-1 

1 a- 

j“0 
m-1 

■ < - r- N b j > a 


- i- V .1 
k o 


-kl Bi ’ 


j-o 

-1 1 i 

As would be expected B --r-y b 

K. aw 

0 j =0 


We will now close this section by proving two propositions which make 
clear the method of solution we have adopted. 

Let MW be the vector space of mxn matrices over the field F. Let 
^n 1,0 tll ° vector space of mnxl vectors over F. Then we have r.ho obvious 


-12 


vector space isomc, iiiwn ft 


11 

>12 


In 


ft 


21 

? 22 


2n 



MW ** M defined ast 



n 


P 11 

P 12 

* ‘ * P ln 

P 21 

P 22 

• 

• 

* * * P 2n 

p ml 

• 

• 

P m2 

* * * P mn 





Let G be as in (3.1). Let the polynomial tt(u) in F(u) be the characteristic 
9 

polynomial of G , Tt(u) ■ det ( u I -G ) . Let II ■ ( ti(u) ) bo the principal 
g mn g 

ideal in F(u]. Then F(u]/II is a ring. 

Define the function ht F(u]/Il— F(x,y)/4' in the following manners 


ht II + a(u)f~> Y + a(g(x,y)) 


Proposition 3.1 . The function h is a ring homomorphism. 

Proof : 

We first show that h is well-def i.ied. I«t II + a(u) - II + b(u) (i.e., 

a(u)-b(u) ■ k(u) tt(u) ) . Show that Y + a(g(x,y)) “ Y + b(g(x,y)) i.c., 

show that a (g (x,y) ) - b(g(x,y') - ^ (x,y) $ 2 (x) + c 2 (x,y) ^ (y) . I claim 

that tt ( g (x,y) ) *» k (x,y) <P 2 (x) tj> 2 (y) + ^ (x,y) <J> 2 (x) + k 2 (x,y) i^ 2 (y) . 

Suppose that we are working over K[x,y] . We have Tt(u' ■ (u-v^^(u-v^ 2 ) . 

(u-v ) where v.,, v. v are the mn eignevalues of G . Let 
mn 11 12 mn g 
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V l j " 9(' A # li j) • Then 

n m 

n<g(x,y)) ■tt n (<! (x,y) 0*3) 

i-1 j-1 

Now wo can show tliat each factor g(x,y) - g(A^#Mj) con he written in the 
form 

g(x,y) - g(X A Uj) - (x,y) (x-A^ ♦ (y) (y-y^) 

This can be seen easily from the fact that If g(x,y) - g fc x t + g t _jX* 1 +••• 
gjX + g Q then 

g(x,y)-* ,-W^ - l9 fc x t-1 ♦ (g t-1 + 9 t A >r t_2 + (g t _ 2 ♦ g^A, + 9 t iJ ,xt 3 ♦ • •• 

+ (9j + 9 2 X 2 + “* 9 t X i + gtAj^.y) - g(A i *P^> 

Therefore (3.3) can be written as 

n m 

Tl(g (x,y) ) -71 TT (kj.j (x*y) (x-A^) + (y) (y-y..)) 

i-1 j-1 

In expanding this product we see that every term in the sum will be 
of either of the two forms, a(x,y)4> 2 (x) or b(x,y)l|/ 2 (y) . Therefore 
Mg(x,y)) - t 1 (x,y)4» 2 (x)<j< 2 (y) + (x, y) <J> 2 (x) + q 2 (x,y)il> 2 (y) in form (2.1) 

over K[x,y]. Since F[x,y] C K(x,y] and form (2.1) is unique we must have 
that tjfxjy), PjtXjy) q^lx^) are actually elements of F(x,y], 
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Thcrefore v (g(x,y) e V and h is woll defined. Now h is a ring homo- 
morphism since 


and 


h[ll + a(u) + II + b(u)] ■ h[n + (a(u) + b(u))] 

* h(Il + a(u)) + h(H ♦ b(u)) 

h [ ( II + a(u) ) ( II + b(u) )] - h( II + a(u)b(u) ) 

■ hJH + a(u))h(.JI + b(u)) 


and h(ft+ 1) - V ♦ 1. 

This completes the proof of Proposition 3.1. 

Nov; since MW i.s an f’tx,y]/H' -module, M a FluI/H-module and h:Ftu]/II ■ > F[x,y]/H' 

n 

a ring homomorphism, MW can be made into an F[u]/H-module in the natural way. 
Define multiplication (*)sF[uJ/lIx MW -► MW by: 


(11+ a(u))*P “ h(n+ a(u))*P 
We now have 

Proposition 3.2 . The map f is an F lul/II -modulo isomorphism- 
Proof : 

We already know that f is a vector space isomorphism. In order to 
show that f is an F [u J/11-module isomorphism we just need to show that 

f (( IT + a(u))*j>) = (II + a(u))*f(£) - (11+ a(u))*P * h ( n + a(u))*P 

Let us shew that 

f((II + u)*jj) - h(It + u)*P. 

Well 
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f«n + u)*p>- f(c g *p) 

- X n J *p*A k 

jk 3 

- (Y ♦ g(x,y) ) *P 
■ h (II + u)*P 

Now it is clear by induction that 

f((Il + u 1 )*^) - h((n + u*))*P 

Therefore 


f ( ( II + a(u))*£) ■ h ( II + a(u))*P 

This completes the proof of Proposition 3.2. ^ 

The equation we are considering takes the two forms given in (3.2) 
and (1.1). We know that in order to obtain the unique solution to (3.2) 
we have to invert the matrix G g or equivalently find the inverse of + u 
in FluJ/ff. The above Propositions show this is the same as obtaining the 
inverse of f + g(x,y) in F[x,y]/Y while working with form (1.1) of the 
equation. 

In the following two sections we will be concerned with the problem 
of constructing the solution to several special cases of the general equa- 
tion. It is of course assumed that a unique solution does exist. We also 
prove a stability theorem associated with the Lyapunov equation. 


4. The Equation PA + BP = Q 

As shown when proving the Main Theorem the solution to equation 


PA + BP = Q is given by 
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P - f (f (x,y)«>od'f »Q) 

nA 

whoro f(x,y) e F(x,yl , (Y + f(x,y))(¥ +(x + y) ) - Y + 1. 

It has also been mentioned that such an f(x,y) can be found by using 
Resultant Theory [ G) or from a constructive proof of the Hilbert-Nullstel- 
lcnsatz. But in simple cases like this wo need not resort to such general 
theory. 

In carrying out computations, it may be advantageous instead of finding 
f(x,y) such that f(x,y)g(x,y) * k^fx) + k^fy) +1 to find f u (x,y) such 
that f u (x,y)g(x,y) * k^tj^fx) + ^ 2 ^ 2 ^ + u w ^ ore u i- s any non-zero element 
in P. The solution P is then given by P=(l/u)*f (f (x,y)mod¥,Q) . 

We construct f^tx.y) in this manner. 

We do have that 


where 


Let 


x + y|<J>^(x)iJ> 2 (y) - 4> x (y) (x) 


Ol<x) * <J> 2 (-x), ^(x) = * 2 (-x). 


p(x,y) 


4> 2 (x)if> 2 (y) - ^(yWjU) 

x+y 


( 4 . 1 ) 


Since <$> 2 (x),^ (x) are coprime jW o for all i,j) we have 

* o (x), P e (x), A. (x) , y^ (x ) 


X e (x)*l<x) + (x) 4> 2 (x) - e 

X^(x)i|( 2 ( x ) + pMxl^lx) 


( 4 . 2 ) 


= e 


- 17 - 


Let f (x,y) - X (x)p ’ (y)p(x,y) , since 
u c c 

f (x,y) (xfy) - x (x)p' ( y)p(x,y) (x+y) 

U GO 

- > (x)M ' (y) <P_ (x) Ip (y) + eX'(y)ip (y) 

+ ep c (x)<{» 9 (x) - U o (x)X^(y)<P 2 (x)<p 2 (y) - e 2 

With u - -e 2 wo have (V + f u (x,y))(Y + (x+y)) ■ ¥ + u. 

A different method for obtaining an f(x,y) such that 
f(x,y)(x+y) * k 1 <P 2 (x) + 2 <y) + 1 is the following: 

Divide <p 2 (x) b V x+ y in x * 

4> 2 (x) ■ h(x,y) (x+y) + h(y). 

For x » -y we have ^(-y) “ 4^ (y) ■ h (y ) . Now since $ (y) ,ip 2 (y) are 
relatively prime there exist X(y),y(y) such that 


^ <y) (y> +M(y)^ 2 (y) ■ i 

^ ^(y)[<P 2 (x) - h(x,y) (x+y) ) +M(y)ip 2 (y) - 1 



-X 


(y)h(x,y) (x+y) + M <y) <J > 2 (x) + U(y)>P 2 (y) 


1 


Let f(x,y) = -)i (y)h(x,y) . 


Tho Lyapunov Equation PA + A'P = Q 

The Lyapunov equation is a special case of PA + BP = Q, B = A', A is 
stable. With the appropriate modifications to the first procedure for 
constructing the solution we have: 
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Lct 1 (x) (x) be such that 
o e 

T (x)<f>. (x) + M (X)4> (x) - e e / 0 

el e i 

4>_ <x) <f> <y) -4> <x) <f» Cy) 

P(x,y) - (4.3) 

f (x,y) » T (x)T (y)p(x,y) 

U C C 


Algorithm for solving the Lyapunov equation A'P + PA ■ Q. 

R^) Obtain <J> 2 (x) the characteristic polynomial of A. 

<|> 2 (x)4> (y) - <f> (y) (x) 

R 2 ) Sot p(x,y) - — . 

R^) Using the Extended Euclidean algorithm or an equivalent 

method obtain T (x) and e. 
e 

R ) Form f (x,y) - T (x)T (y)P(x,y). 

4 u e e 

R 5 ) Find P u - f BA (f u (x,y)mod4',Q) . 

1 2 
R,. ) Set P = — P , u = -e . 

6 u u 

Computer Implementation 

Since we are interested in an exact computer solution we restrict 
the field of interest to that of the rational numbers R. The algorithm 
is fully implementablo, using the remarkable facilities provided by the 
computer programming system MACSYMA available at M.I.T. MACSYMA is a 
large computer programming system used for performing symbolic as well 


as numerical computations. 
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Three versions of the algorithm have been constructed and programmed 
on MACSYMA. They are the Rational Algorithm, the Integer Algoritlim, and 
the Modular Algorithm having names indicative of the mode in which arith- 
metic operations are carried out. 

The Rational Algorithm 

It consists of carrying out steps R, through R, in rational arithmetic. 

L D 

The Integer Algoritlun 

Suppose that matrices A and Q contained integer entries. The poly- 
nomials $ 2 (x), p(x,y) then have integer coefficients. 

. . , . . n n-1 

Let d)-(x) ■ a x + a ,x ■ , . . + a_ 

< n n-i 0 

< ^l ( x ) " V" + Vl^' 1 + + d 0 

Define the 2nx2n matrix S 



a 

0 

0 

d 

0 

0 


n 



n 




a , 

a 

0 

d , 

d , 

0 


n-1 

n 


n-1 

n-1 



• 

a , 

• 





• 

n-1 

• 





• 


• 






• • • 

a 

d. 

d 

d 




n 

1 

2 

n 

s » 

a„ 

a. 

a , 

d. 

• • • 

d , 


0 

1 

n-1 

0 


n-1 


0 

a_ 

a 

0 

d 




0 

n-2 


0 



0 

0 

• 

0 

0 

• 

• 




• 



• 




• 

• 

• 






• 

• 



• 

• 


• 

• 



• 

• 






• 

• 






0 

0 


0 

0 

d 




0 



0 
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Wo know that dots f 0 since it is the resultant of ^(x) and $ 2 (x) 
which are coprime. If we lot o * dots the linear system 


r x .1 

n-1 


”o " 

X n-2 


0 

X 0 



T n-1 

- 


T n-2 



_ T 0 


_e _ 


has an integer solution and we have integer polynomials 

T (x) ■ T . x n 1 + ... T_, A (x) “ A ,x n * + ... + A„ which satisfy 
e n-1 0 e n-1 0 J 

T 0 ( x )f(x) + A (x) » e 

e l e 2 

This means that in (4.3) has integer coefficients and so does 

f (x.yjmod'i', which implies that P ** f (f (x,y)mod , J',Q) has integer 

U \2 BA VI 

entries. 

The algorithm proceeds as follows. 

1^) Find ^(x) t * ie characteristic polynomial of A. 

4-2< x) <J> 2 (y) ” ♦l (x) *l (y) 

I 2 ) Set p(x,y) = — . 

I,) Find t (x) and e. 

J c 
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I ) Form f (x,y) - t (x)i (y)p(x.y). 
4 U CO 

I ) Find P » f (f (x, y)modY,Q) . 

D U dA U 

I,) Set P ■ — *P , u “ -e^. 

6 u u 


The Modular Algorithm 

The integer algorithm paves the way for a modular approach to the 
solution. 

Suppose p is a prime that does not divide e « detS. If A ■ (a^) 

and Q “ (q^) aro matrices with integer entries let p Q » (a^modp) and 

A “ (a. .modp) be considered as matrices over Z . A left subscript p 
P ij P 

on a polynomial b(x,y) written as p b(x,y) denotes coefficient reduction 
modulo p. Suppose that coefficient arithmetic is done modulo p. We then 
have 


p <P 2 (x) - detdx - p A) 


J?(x,y) 


^ 2 <x) <D 2 (y) - ^(x) <*» 1 (y) 


x+y 


pWi 1 * 1 ♦ 


p X e (x) p * 2 (x) 


e 

P 


f 

P u 


(x,y) 


T (x) T 
p e p e 


(y) p P(x,y) . 


Let P = f_„ (f (x,y)mod V, Q) where all arithmetic is done modulo p and 
p u BA u 1 p p x * 

V = ( d>_ (x) , <{>„(y)) in Z (x,y]. If P and u are obtained for a suf- 

P P 2 P 2 p pup 

ficient number of primes, the Chinese Remainder Theorem can be used to 

find P and u making it possible to obtain the solution P = — *P . 
u u u 
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Tho algorithm is an follows: 


V 

Obtain A, Q. 

p p* 


M 2> 

Obtain d>_(x) ■ dct(Ix - A). 

P 2 p 


M 3> 

. „V X W' ,) - „V x W vl 

Set P(x,y) - *■ c *- . 

p ' x+y 

• 

V 

Obtain T (x) , e. 
p e p 


V 

• V 

Set f (x,y) ■ T (x) T (y) P( x »y)» 
p u p e p e p 


V 

Obtain P ■ f„„( f (x,y)mod Y, Q) . 
p u BA p u 1 p p x 


M 7> 

Repeat steps 1-6 for a sufficient number of 

primes and 


using the Chinese Remainder Theorem find 

and u “ -e 

V 

Set P “ -»P . 

u u 



Since considerable coefficient growth takes place in intermediate 
computations of tne Integer Algorithm, a lot of storage is being used 
up. In such cases it is advantageous to use the Modular Algorithm. 

Arithmetic Complexity of the Integer Algorithm 

We are concerned with the number of integer operations (addition, 
subtraction, multiplication, division) involved in running the Integer 
Algorithm when A and Q are nxn matrices, using classical operations. 


Step 1^ : There are several methods for obtaining the Characteristic 

polynomial <(> 2 (x) of a stable matrix. Evaluating 4> 2 (x) at n distinct 
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4 

points and then solving for thu coefficients requires 0(n ) operations. 

If n is small (say n < 20), evaluating ({^(x) at x ■ 1 where $ 2 (1) - A, 

X “ (log^Al and then at x « 10* allows ono to “read off" the coefficients 
of f{> 2 (x) from a large integor. This procedure requires only 0(n 3 ) 
operations. 

2 

Step I ^ '• This step can be done in 0(n‘) operations. 

Step I^i Solving a linear set of 2n equations simultaneously is and 

3 

0 (n ) operation. 

Step It Performing the multiplication as I (x) (I (y)*P(x,y)J requires 

*i © 

3 

0(n ) operations. 

Step I^r Obtaining f^ (x,y)mod4' involves two polynomial divisions 
can be done in 0(n 3 ) operations. To form f (f (x,y)modY,Q) we use 

BA U 

4 

0(n ) operations. In the event that the matrix Q is a product of vec- 
tors Q - c»c* this calculation can be done in 0(n 3 ) operations. 

2 

Step I,: It can be done in 0(n ) operations. 

o 

It can therefore be seen that the overall calculation requires 
4 , 3 

0(n ) operations in general aid 0(n ) operations in the special cases 
mentioned. 

Storage requirements are much harder to determine since the imple- 


mentation is on a variable length word computer. 
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Numerical Examples 

We now continue this section by gi\ .ng two numerical examples . 
We wish to compute 

QO 

G - J x* (t) *Q. x(t)dt 
where x(t) is a solution to 

x (t) - Ax(t) x (0) - c (*) 

.Jr 

The system modelled by (*) is of the form 



where the number of blocks is finite. 

Example It The number of blocks is 5 with £*1, K**l, M=10000. Listed 
are the corresponding A matrix, the Q matrix and the solution P to the 
equation PA + A'P = Q. 

Example 2: The number of blocks is 2. Listed are the corresponding 

A matrix in parametric form, the Q matrix and the parametric solution 
P of the equation PA + A'P = Q. The parametric solution P is valid 
only for appropriate values of E, M, Z, (Z=0 . 


///////////, 


25 - 


0 


0 


0 


0 


0 


26 


_1 

7 

0 

0 

0 

0 

0 

_12500 

3 

0 

_10000 

3 

0 

0 

0 

_1 

2 

0 

0 

0 

_ 10000 
3 

4 

0 

20000 

3 

0 

0 

0 

0 

0 

_1 

2 

0 

-2500 

0 

-5000 

0 

0 

0 

0 

0 

0 

A 

. _5000 

A 

_10000 

n 

u 

3 

U 

3 

W 

0 

0 

0 

0 

0 

A 

_2500 

ft 

_ 5000 

ft 

0 

3 

u 

3 

1/ 


0 

0 

0 

0 

0 

-2500 

0 

_5000 

3 

0 

_2500 

3 

0 

• 0 

0 

0 

0 

-5000 

ft 

10000 

ft 

5000 

U 

» 

U 

3 

0 

0 

0 

0 

0 

-7500 

0 

-5000 

0 

-2500 

0 

„1 

1 

0 

0 

0 

-5000 

0 

_20000 

3 

0 

_10000 

3 

0 

0 

Or 

_1 

~2 

0 

-2500 

0 

_10000 

3 

0 

_12500 

3 


Exampla 2. 


0 10 0 

2E _ 2Z E Z 

M H M M 

0 0 0 1 

E , Z _2E _2Z 

M M M M 



0 0 
0 0 
0 0 
0 1 


P ** 


K 

2Z 

0 

0 


0 

3Z 

0 

6Z 


- E 

2Z 


0 

_M_ 

6Z 

0 

_M_ 

3Z 


0 


0 
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In closing this section it is interesting to observe how the ideas 
presented here can be used to prove stability theorems constructs /ely 
In particular we prove 

Theorem 4.1 . Let A bo an nxn matrix over the reels R. Let C be a 
pxn matrix. If A is a stability matrix and (A,C) an observable pair 
then the equation PA ♦ H'P • -C'C has a unique symmetric positive definite 
solution P. 

Before proceeding with the proof we introduce the notion of a positive 
polynomial in R(x,y] [2], Jfp(x,y) is in R(x,y] we can write it 
os 

p(x,y) - f.' (y)C(p)Mx) 

where f.(z) is the column vector l,z,...z n 1 , n with n being one plus 

the largest power of p(x,y) in either x or y and C(p) an nxn matrix 

over R. This introduces a bijection between R{x,y] and the set of all 

square matrices. Wo then call a polynomial p(x,y) positive if C(p) is 

i) symmetric and ii) positive definite. One can then prove l 2] that 

p(x,y) in F x,y] is positive if and only if thero exist polynomials 

it, * • •# * (m the size of C(p)) such that 

12 m 

m 

. . V* n (x)7J (y) 

p(x,y) - 1 i 

i-1 

where {n^ (x> } are a basis for the vector space (over R) of polynomials 
of degree less than m. One can also prove I 2 I that the f u (x,y) mod ¥ 
given in (4.3) is positive. 
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Pr oof of Tl» or«‘W 4 . 1 

Since A Is a stability matrix A^ ♦ A^ / 0 for all i,j therefore 
a unique solution P to tho equation I . ♦ A'P ■ -C'C exists. Wt? also have 
that f u (x,y) mod ¥ is positive. We can therefore write 

f u <x,y) mod ¥ « 

Ke know that the unique solution is given by 


£V«>w l( y) 


P - — f (f (x,y) mod¥, -C'C) (u » -e ) 

u DA u 

‘ V w£V’‘ , V'' ) - C ' C > 

e 1=1 

- — ^tt^A'JC'C tt.(a) 

e i“-i 

Since C'C > 0 we have P 0. 

Suppose now that there exists z = 0 such that z'Pz=0. 



1 |CTt. (A) z I I = 0 

for 

l<i<n 


cr i ( a) z = o 

for 

l<i<n 


Since {tk} are a basis there exists a matrix T such that 


’tt (x)' 


1 

1 


X 

TT- (*) 



2 

rs 

• 

• 


• 

TT (X) 


n-1 

X 

L n J 


™ — 
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i-1 


> f BA (t il' 7 , l (x) + ,,,+ t in’ lf n (x) ' C) " CA l l i .l n whcre (t n' * ' * * t i n ) 
the i^h r ow of T 


> Jyt ctoo " ca11 




Oefine the operator H » R n R nr> by: 


C 

CA 


H (w) = 


CA 


n -1 


w 


Since (A,C) observable pair the null space of H is {o}, and since 

n 

i (A)z ■ 0 this implies that J^t . . 7T . (A) z - 0 for all l<i<»\ 

j-1 13 3 


Cv 


-> H(z) - 0. 

This is a contradiction since z^O. 


5 . The Equation P-BPA = Q 

We again wi3h to construct f(x,y) such that (4* + f(x,y)KV + (1-xy) ) «* 
V+l. Let 

<t> (x) = det(Ix-A) =• a x n + a .x n 1 +...+ a,x + a 
2 n n-J. i u 

il/_ (x) «* det(Ix-B) = b x m + b .x m ^ +. . . b x + b 
Y 2 m m-1 1 0 

. . n n-1 

(x) = a Q x + a^x +...+ a^ 

V X> " V" + b l X ' n ' 1 + -** + b m 

From the above definition we can see that the roots of <J> 3 (x) are the 

values -p — where X. / 0. Since we assume that a unique solution exists 
A i 1 
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we must have 1- A^ / 0 r ot all i,j. Then we must have that $ (x) , 
ij'^(x) are coprimo. Because if they have a non-trivlul factor k(x) they 
must also have at least one common root (i.e. # A^ - for at least 
some (i,j)). 

On the other hand we also have that 


a) if n >m then 1 - xy | y n_m «^ 2 (x) ^ (y) - <> 3 (y) ip 3 (x) 

b) if n < m then l-xy|x ni n <|> 2 (y) <J/ 2 (y) - <t> 3 (y) tp 3 (x) 

- ^ 

We are now ready to construct f(x,y). 

Since 4> 3 (x) ^ (x) are relatively prime wo have A(x) p(x) A' (x) A'(x) 
such that 


A (x) i^ 3 (x) + p (x) <t> 2 (x) « 1 


A • ( x ) ( x ) + p*(x)<J> 3 (x) = i 


If n>m let 


y n " m 4> 2 (x)^ 2 (y) - <j> 3 (y)ip 3 (x) 
1 - xy 


p(x f y) 


if n<m let 


m-n , 


p(x,y) 


x <P 2 <x) 4> 2 (y) - <f> 3 (y)ij> 3 (x) 
1 - xy 


Then f(x,y) = p(x)p'(y) p(x,y). 

The discrete Lyapunov equation P-A'PA = Q is a special case. 

6. Over Integral Domains 

Suppose now that we are investigating equation (1.1) over E some 
integral domain. The next proposition gives a necessary and sufficient 
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condition for the existence of a unique solution to (1.1) for all Q. 

Proposition 6.1 . Equation (1.1) has a unique solution over E, for each 
Q, iff H* + g(x,y) is a unit in Elx^l/T. 


Proof ; 

Let P ■ f BA (f ( x »y) modl’.Q) where f(x,y)g(x,y) ■ (x,y) <J> 2 (x) + 

k 2 (x,y) (y) 41 







f BA ,g(X ' yJ ' P > 
f BA (g(X ' y) ' f BA (f(X,y) ' Q)) 


W 1 ' 2 ) 


The solution P is unique. This follows in the same manner as in the 
proof of the Main Theorem. 

Suppose that equation (1.1) does have a unique solution for all Q. 
This means that G in (3.2) is invertible. We have that Tr(u) = det(Iu-G ) 

g g 

From the Caylcy-Hamilton theorem if tt (u) =* tt^u*' + tt^ ^u*" 1 +...+ it 


TT(G J = TT G +...+ it I 
g t g 0 


T . C i \ t t-1 

Let f (u) = - — u 
IT. 


TT 


t-1 t-2 

u 

TT, 


- — . Then f (u) • u + — . IT (u) =1 
TT TT 

0 0 


0 0 

Therefore II + u is a unit in E[u]/Il. The proof of Proposition 3.1 


remains valid. Therefore 


(11+ f (u)) (n + u) = n + 1 
* h(*+ f (u)) • h(n + u) = T +1 

-> (Y + f (g (x,y) ) (V + g(x,y)) = V + 1 

which means that V + g(x,y) is a unit in Elx^l/f. 
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